#*************        丁香园-NHANES      *************
#*************        Analyze Code       *************
#*************       Table2 专项讲解     *************
#*************      未经允许不准转载     *************

#### 0.读取数据 & 准备环境 ####
library(gtsummary)
library(survey)
library(tidyverse)
setwd("F:/Rproject/课程_5.4_Table2_专项讲解以及演示")

paper.data <- readRDS('Table2_Sample_Data.RData') # 读取课后材料中准备的数据

# tab3 有Quartile of total L and Z 计算
paper.data$Total.LZ.quantile.var <- cut(paper.data$Total.LZ,
                                          breaks = quantile(paper.data$Total.LZ),
                                          labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
# 加权情况相下测试
NHANES_design <- svydesign(data = paper.data, ids = ~SDMVPSU, strata = ~SDMVSTRA, 
                           nest = TRUE, weights = ~WTDRD1, survey.lonely.psu = "adjust") 

Total.LZ.quantile.res <- svyquantile(~ Total.LZ, NHANES_design, quantiles = c(0, 0.25, 0.5, 0.75, 1))

paper.data$Total.LZ.quantile.var <- cut(paper.data$Total.LZ,
                                          breaks = Total.LZ.quantile.res$Total.LZ[,'quantile'],
                                          labels = c('Q1', 'Q2', 'Q3', 'Q4'))

NHANES_design <- svydesign(data = paper.data, ids = ~SDMVPSU, strata = ~SDMVSTRA, 
                           nest = TRUE, weights = ~WTDRD1, survey.lonely.psu = "adjust") 

#### 1.Table2 基础版本，获得表格主体 ####
##### 1.1 支持的回归模型列表 ##### 
# stats::lm(), stats::glm(), stats::aov(), ordinal::clm(), ordinal::clmm(), survival::coxph(), 
# survival::survreg(), survival::clogit(), lme4::lmer(), lme4::glmer(), lme4::glmer.nb(), 
# brms::brm(), geepack::geeglm(), gam::gam(), mgcv::gam(), nnet::multinom(), survey::svyglm(), 
# survey::svycoxph(), survey::svyolr(), MASS::polr(), MASS::glm.nb(), mice::mira, lavaan::lavaan(),
# cmprsk::crr(), stats::nls(), lfe::felm(), rstanarm::stan_glm(), VGAM::vglm(), and more!

##### 1.2 参数说明 ##### 
# exponentiate: Logical indicating whether to exponentiate the coefficient estimates. Default is FALSE.
# 上面的那个参数是发文以指数形式呈现 
# include: Table 中要呈现的变量，以及变量的排序
# label: 手动更改 Table 中变量的呈现内容，例如将 Age 这个变量在 Table 中呈现为 Age (years)

# conf.int: 是否显示置信区间，Logical indicating whether or not to include a confidence interval in the output. Defaults to TRUE.
# conf.level: 置信区间范围，默认为 0.95

# 使用函数 tbl_regression，来一键生成结果
# https://www.danieldsjoberg.com/gtsummary/articles/tbl_regression.html
# Y1 = CERAD.delay.recall
# Fully adjusted

  

m1 <- svyglm(CERAD.delay.recall ~ Dietary.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
             design = NHANES_design)
tbl_regression(m1)
# 删除p值列 合并部分列值显示
aaa <-tbl_regression(m1, exponentiate = TRUE, include = c(Dietary.LZ))%>%# 只呈现 Dietary.LZ
  modify_column_hide(p.value)%>% # 隐藏p值
  modify_header(estimate ~ 'exp(Beta) (95% CI)')%>% # 编辑表头
  modify_column_merge(pattern = "{estimate} ({ci})")  # 合并列值
# 使用之前先看下 show_header_names进行观测数据列进行合并
tbl_regression(m1, 
               exponentiate = TRUE, 
               include = c(Dietary.LZ, Age, Sex),
               label = list(Age ~ 'Age (years)', Sex ~ 'Sex'),
               conf.level = 0.80) # 呈现多个变量 Dietary.LZ


#### 2.Table2 进阶版本（一）：加上列：模型统计量等信息 ####
# add_glance_table(): 加入模型的统计量，包括 r.squared, AIC, sigma 等
# 可以在默认输出的统计量中，选择自己所需要的，不同的模型得到的统计量不同，详见 help
# add_global_p(): 分类变量的汇总 P value
# inline_text.tbl_regression()：以特定的模式获取表格中的模型结果，
# 默认模式为："{estimate} ({conf.level*100}% CI {conf.low}, {conf.high}; {p.value})"

tbl_regression(m1, 
               # exponentiate = TRUE, 
               include = c(Dietary.LZ, Age, Sex, Alq.group),
               label = list(Age ~ 'Age (years)', Sex ~ 'Sex'),
               conf.level = 0.95,
               pvalue_fun = function(x) style_pvalue(x, digits = 2)) %>% 
  add_global_p() %>% # 加上分类变量的显著性 p value
  modify_header(
    # 这样可以添加列 我以为是add_xx 发现不行 2023年3月25日19:10:20
    statistic = "**Statistic**",
    std.error = "**SE**"
  )%>%
  add_glance_table(include = c(AIC, nobs)) %>% # 加上模型的统计值
  modify_footnote(ci = "CI = 修改模型的脚注 CI", abbreviation = TRUE) %>% # 修改脚注
  bold_labels() %>% # 加粗 label 的结果
  italicize_levels() %>% # 斜体化 level
  bold_p(0.05)  # 加粗 p 显著的结果
  

#### 3.Table2 进阶版本（二）：左右叠加不同的 Table ####
# Dietary.LZ
dr.age.adjusted.tbl <- svyglm(CERAD.delay.recall ~ Dietary.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Dietary.LZ),
  label = list(Dietary.LZ ~ 'Age-adjusted'))



dr.fully.adjusted.tbl <- svyglm(CERAD.delay.recall ~ Dietary.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                          design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'Fully-adjusted'))

# Supplement.LZ
supp.age.adjusted.tbl <- svyglm(CERAD.delay.recall ~ Supplement.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Supplement.LZ),
                 label = list(Supplement.LZ ~ 'Age-adjusted'))


supp.fully.adjusted.tbl <- svyglm(CERAD.delay.recall ~ Supplement.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                          design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Supplement.LZ),
                 label = list(Supplement.LZ ~ 'Fully-adjusted'))

# Total L.Z gcr
Total.age.adjusted.tbl <- svyglm(CERAD.delay.recall ~ Total.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'Age-adjusted'))


Total.fully.adjusted.tbl <- svyglm(CERAD.delay.recall ~ Total.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                  design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'Fully-adjusted'))

# Total L.Z Quartile gcr

Totallz.age.Quartile.adjusted.tbl <- svyglm(CERAD.delay.recall ~ Total.LZ.quantile.var + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'Age-adjusted'))


Totallz.fully.adjusted.tbl <- svyglm(CERAD.delay.recall ~ Total.LZ.quantile.var + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                   design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'Fully-adjusted'))


# regression 的结果一般采用纵向拼接（stack）而不是横向拼接（merger）
CERAD.delay.recall<-tbl_stack(
  tbls = list(dr.age.adjusted.tbl, dr.fully.adjusted.tbl,
              supp.age.adjusted.tbl, supp.fully.adjusted.tbl,Total.age.adjusted.tbl,Total.fully.adjusted.tbl,Totallz.age.Quartile.adjusted.tbl,Totallz.fully.adjusted.tbl),
  group_header = c("Dietary.LZ", "Dietary.LZ", "Supplement.LZ", "Supplement.LZ","Total.LZ", "Total.LZ","Quartile of total L and Z","Quartile of total L and Z")) 



# CERAD1

dr.age.adjusted.tbl <- svyglm(CERAD1 ~ Dietary.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'Age-adjusted'))


dr.fully.adjusted.tbl <- svyglm(CERAD1 ~ Dietary.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'Fully-adjusted'))

# Supplement.LZ
supp.age.adjusted.tbl <- svyglm(CERAD1 ~ Supplement.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Supplement.LZ),
                 label = list(Supplement.LZ ~ 'Age-adjusted'))


supp.fully.adjusted.tbl <- svyglm(CERAD1 ~ Supplement.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                  design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Supplement.LZ),
                 label = list(Supplement.LZ ~ 'Fully-adjusted'))

# Total L.Z gcr

Total.age.adjusted.tbl <- svyglm(CERAD1 ~ Total.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'Age-adjusted'))


Total.fully.adjusted.tbl <- svyglm(CERAD1 ~ Total.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                   design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'Fully-adjusted'))

# Total L.Z Quartile gcr

Totallz1.age.Quartile.adjusted.tbl <- svyglm(CERAD1 ~ Total.LZ.quantile.var + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'Age-adjusted'))



Totallz1.fully.adjusted.tbl <- svyglm(CERAD1 ~ Total.LZ.quantile.var + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                     design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'Fully-adjusted'))

# regression 的结果一般采用纵向拼接（stack）而不是横向拼接（merger）
CERAD1col<-tbl_stack(
  tbls = list(dr.age.adjusted.tbl, dr.fully.adjusted.tbl,
              supp.age.adjusted.tbl, supp.fully.adjusted.tbl,Total.age.adjusted.tbl,Total.fully.adjusted.tbl,Totallz1.age.Quartile.adjusted.tbl,Totallz1.fully.adjusted.tbl),
  group_header = c("Dietary.LZ", "Dietary.LZ", "Supplement.LZ", "Supplement.LZ","Total.LZ", "Total.LZ","Quartile of total L and Z","Quartile of total L and Z")) 



# cerad2
dr2.age.adjusted.tbl <- svyglm(CERAD2 ~ Dietary.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'Age-adjusted'))

dr2.fully.adjusted.tbl <- svyglm(CERAD2 ~ Dietary.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'Fully-adjusted'))
# Supplement.LZ
supp2.age.adjusted.tbl <- svyglm(CERAD2 ~ Supplement.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Supplement.LZ),
                 label = list(Supplement.LZ ~ 'Age-adjusted'))

supp2.fully.adjusted.tbl <- svyglm(CERAD2 ~ Supplement.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                  design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Supplement.LZ),
                 label = list(Supplement.LZ ~ 'Fully-adjusted'))
# Total L.Z gcr
Total2.age.adjusted.tbl <- svyglm(CERAD2 ~ Total.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'Age-adjusted'))

Total2.fully.adjusted.tbl <- svyglm(CERAD2 ~ Total.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                   design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'Fully-adjusted'))

# Total L.Z Quartile gcr
Totallz2.age.Quartile.adjusted.tbl <- svyglm(CERAD2 ~ Total.LZ.quantile.var + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'Age-adjusted'))

Totallz2.fully.adjusted.tbl <- svyglm(CERAD2 ~ Total.LZ.quantile.var + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                      design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'Fully-adjusted'))
# regression 的结果一般采用纵向拼接（stack）而不是横向拼接（merger）
CERAD2col<-tbl_stack(
  tbls = list(dr2.age.adjusted.tbl, dr2.fully.adjusted.tbl,
              supp2.age.adjusted.tbl, supp2.fully.adjusted.tbl,Total2.age.adjusted.tbl,Total2.fully.adjusted.tbl,Totallz2.age.Quartile.adjusted.tbl,Totallz2.fully.adjusted.tbl),
  group_header = c("Dietary.LZ", "Dietary.LZ", "Supplement.LZ", "Supplement.LZ","Total.LZ", "Total.LZ","Quartile of total L and Z","Quartile of total L and Z")) 

# CERAD 3

dr3.age.adjusted.tbl <- svyglm(CERAD3 ~ Dietary.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'Age-adjusted'))
dr3.fully.adjusted.tbl <- svyglm(CERAD3 ~ Dietary.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'Fully-adjusted'))
# Supplement.LZ
supp3.age.adjusted.tbl <- svyglm(CERAD3 ~ Supplement.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Supplement.LZ),
                 label = list(Supplement.LZ ~ 'Age-adjusted'))

supp3.fully.adjusted.tbl <- svyglm(CERAD3 ~ Supplement.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                  design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Supplement.LZ),
                 label = list(Supplement.LZ ~ 'Fully-adjusted'))
# Total L.Z gcr
Total3.age.adjusted.tbl <- svyglm(CERAD3 ~ Total.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'Age-adjusted'))

Total3.fully.adjusted.tbl <- svyglm(CERAD3 ~ Total.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                   design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'Fully-adjusted'))
# Total L.Z Quartile gcr
Totallz3.age.Quartile.adjusted.tbl <- svyglm(CERAD3 ~ Total.LZ.quantile.var + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'Age-adjusted'))
Totallz3.fully.adjusted.tbl <- svyglm(CERAD3 ~ Total.LZ.quantile.var + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                      design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'Fully-adjusted'))

# regression 的结果一般采用纵向拼接（stack）而不是横向拼接（merger）
CERAD3col<-tbl_stack(
  tbls = list(dr3.age.adjusted.tbl, dr3.fully.adjusted.tbl,
              supp3.age.adjusted.tbl, supp3.fully.adjusted.tbl,Total3.age.adjusted.tbl,Total3.fully.adjusted.tbl,Totallz3.age.Quartile.adjusted.tbl,Totallz3.fully.adjusted.tbl),
  group_header = c("Dietary.LZ", "Dietary.LZ", "Supplement.LZ", "Supplement.LZ","Total.LZ", "Total.LZ","Quartile of total L and Z","Quartile of total L and Z")) 



alllei<-tbl_merge(tbls = list(CERAD.delay.recall,CERAD1col,CERAD2col,CERAD3col),tab_spanner = c('CERAD: Score DelayedRecall','CERAD: Trial 1 Score','CERAD: Trial 2 Score','CERAD: Trial 3 Score'))
alllei%>%as_flex_table()%>%flextable::save_as_html(path = "2.8table3.html")

#### Table 3 结束####

#### Table 4 开始####
# Adjusted* beta coefficients (standard error (SE), P-value) for score on CERAD Word Learning sub-test, Animal Fluency
#test, and DSST, for each mg/day increase in L and Z intake, stratified by race/ethnicity

# Dietary.LZ

tb4dr.age.adjusted.tbl <- svyglm(CERAD.total ~ Dietary.LZ + Age, design = NHANES_design)%>%
  tbl_regression(include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'Age-adjusted'),
                 pvalue_fun = function(x) style_pvalue(x, digits = 2)
                 )%>%
  modify_header(
    std.error = "**SE**"
  )
tb4dr.fully.adjusted.tbl <- svyglm(CERAD.total ~ Dietary.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                design = NHANES_design) %>%
  tbl_regression( include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
# Supplement.LZ
tb4supp.age.adjusted.tbl <- svyglm(CERAD.total ~ Supplement.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Supplement.LZ),
                 label = list(Supplement.LZ ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
tb4supp.fully.adjusted.tbl <- svyglm(CERAD.total ~ Supplement.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                  design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Supplement.LZ),
                 label = list(Supplement.LZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
# Total L.Z gcr
tb4Total.age.adjusted.tbl <- svyglm(CERAD.total ~ Total.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
tb4Total.fully.adjusted.tbl <- svyglm(CERAD.total ~ Total.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                   design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
# Total L.Z Quartile gcr
tb4Totallz.age.Quartile.adjusted.tbl <- svyglm(CERAD.total ~ Total.LZ.quantile.var + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
tb4Totallz.fully.adjusted.tbl <- svyglm(CERAD.total ~ Total.LZ.quantile.var + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                     design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'Fully-adjusted'))
# regression 的结果一般采用纵向拼接（stack）而不是横向拼接（merger）
tb4CERAD.total<-tbl_stack(
  tbls = list(tb4dr.age.adjusted.tbl, tb4dr.fully.adjusted.tbl,
              tb4supp.age.adjusted.tbl, tb4supp.fully.adjusted.tbl,tb4Total.age.adjusted.tbl,tb4Total.fully.adjusted.tbl,tb4Totallz.age.Quartile.adjusted.tbl,tb4Totallz.fully.adjusted.tbl),
  group_header = c("Dietary.LZ", "Dietary.LZ", "Supplement.LZ", "Supplement.LZ","Total.LZ", "Total.LZ","Quartile of total L and Z","Quartile of total L and Z")) 

# Animal Fluency score
Animaldr.age.adjusted.tbl <- svyglm(Animal.Fluency ~ Dietary.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
Animaldr.fully.adjusted.tbl <- svyglm(Animal.Fluency ~ Dietary.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
# Supplement.LZ
Animalsupp.age.adjusted.tbl <- svyglm(Animal.Fluency ~ Supplement.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Supplement.LZ),
                 label = list(Supplement.LZ ~ 'Age-adjusted'))
Animalsupp.fully.adjusted.tbl <- svyglm(Animal.Fluency ~ Supplement.LZ + Age + Sex + BMI + Alq.group + 
                                    Smoke.group + PIR + Education.attainment, 
                                  design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Supplement.LZ),
                 label = list(Supplement.LZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )

# Total L.Z gcr
AnimalTotal.age.adjusted.tbl <- svyglm(Animal.Fluency ~ Total.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
AnimalTotal.fully.adjusted.tbl <- svyglm(Animal.Fluency ~ Total.LZ + Age + Sex + BMI + Alq.group + 
                                     Smoke.group + PIR + Education.attainment, 
                                   design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
# Total L.Z Quartile gcr
AnimalTotallz1.age.Quartile.adjusted.tbl <- svyglm(Animal.Fluency ~ Total.LZ.quantile.var + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
AnimalTotallz1.fully.adjusted.tbl <- svyglm(Animal.Fluency ~ Total.LZ.quantile.var + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                      design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
# regression 的结果一般采用纵向拼接（stack）而不是横向拼接（merger）
Animal.Fluencytab<-tbl_stack(
  tbls = list(Animaldr.age.adjusted.tbl, Animaldr.fully.adjusted.tbl,
              Animalsupp.age.adjusted.tbl, Animalsupp.fully.adjusted.tbl,AnimalTotal.age.adjusted.tbl,AnimalTotal.fully.adjusted.tbl,AnimalTotallz1.age.Quartile.adjusted.tbl,AnimalTotallz1.fully.adjusted.tbl),
  group_header = c("Dietary.LZ", "Dietary.LZ", "Supplement.LZ", "Supplement.LZ","Total.LZ", "Total.LZ","Quartile of total L and Z","Quartile of total L and Z")) 


# Digit Symbol Score
#Dietary L and Z
dss.age.adjusted.tbl <- svyglm(DSST ~ Dietary.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )

dss.fully.adjusted.tbl <- svyglm(DSST ~ Dietary.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                 design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )

# Supplement.LZ
dsssupp2.age.adjusted.tbl <- svyglm(DSST ~ Supplement.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Supplement.LZ),
                 label = list(Supplement.LZ ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
dsssupp2.fully.adjusted.tbl <- svyglm(DSST ~ Supplement.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                   design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Supplement.LZ),
                 label = list(Supplement.LZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
# Total L.Z gcr
dssTotal2.age.adjusted.tbl <- svyglm(DSST ~ Total.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
dssTotal2.fully.adjusted.tbl <- svyglm(DSST ~ Total.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                    design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )

# Total L.Z Quartile gcr
dssTotallz2.age.Quartile.adjusted.tbl <- svyglm(DSST ~ Total.LZ.quantile.var + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )

dssTotallz2.fully.adjusted.tbl <- svyglm(DSST ~ Total.LZ.quantile.var + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                      design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
# regression 的结果一般采用纵向拼接（stack）而不是横向拼接（merger）
DSSTcol<-tbl_stack(
  tbls = list(dss.age.adjusted.tbl, dss.fully.adjusted.tbl,
              dsssupp2.age.adjusted.tbl, dsssupp2.fully.adjusted.tbl,dssTotal2.age.adjusted.tbl,dssTotal2.fully.adjusted.tbl,dssTotallz2.age.Quartile.adjusted.tbl,dssTotallz2.fully.adjusted.tbl),
  group_header = c("Dietary.LZ", "Dietary.LZ", "Supplement.LZ", "Supplement.LZ","Total.LZ", "Total.LZ","Quartile of total L and Z","Quartile of total L and Z")) 
# regression 的结果一般采用纵向拼接（stack）而不是横向拼接（merger）

alllei4<-tbl_merge(tbls = list(tb4CERAD.total,Animal.Fluencytab,DSSTcol),tab_spanner = c('CERAD: Total score','Animal Fluency score','Digit Symbol Score'))
alllei4%>%as_flex_table()%>%flextable::save_as_html(path = "2.8table4.html")

#### tab4 结束####


#### tab5 开始####

options( survey.lonely.psu = "adjust" )
 
# Dietary L and Z
CERAtotal.White.adjusted.tbl <- svyglm(CERAD.delay.recall ~ Dietary.LZ+ Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment,design = subset(NHANES_design, Race == 'Non-Hispanic White')) %>%
  tbl_regression(exponentiate = TRUE, include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'White'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")

CERAtotal.black.adjusted.tbl <- svyglm(CERAD.delay.recall ~ Dietary.LZ+ Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment ,design = subset(NHANES_design, Race == 'Non-Hispanic Black')) %>%
  tbl_regression(exponentiate = TRUE, include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'Black'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
# Total L and Z
CERAtotalLZ.White.adjusted.tbl <- svyglm(CERAD.delay.recall ~ Total.LZ+ Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment,design = subset(NHANES_design, Race == 'Non-Hispanic White')) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'White'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
CERAtotalLZ.black.adjusted.tbl <- svyglm(CERAD.delay.recall ~ Total.LZ+ Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment ,design = subset(NHANES_design, Race == 'Non-Hispanic Black')) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'Black'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
#Quartile of total L and Z
CERAtotalQUAN.White.adjusted.tbl <- svyglm(CERAD.delay.recall ~ Total.LZ.quantile.var+ Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment,design = subset(NHANES_design, Race == 'Non-Hispanic White')) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'White'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
CERAtotalQUAN.black.adjusted.tbl <- svyglm(CERAD.delay.recall ~ Total.LZ.quantile.var+ Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment ,design = subset(NHANES_design, Race == 'Non-Hispanic Black')) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'Black'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
table5sdr <- tbl_stack(
  tbls = list(CERAtotal.White.adjusted.tbl, CERAtotal.black.adjusted.tbl,CERAtotalLZ.White.adjusted.tbl,CERAtotalLZ.black.adjusted.tbl,CERAtotalQUAN.White.adjusted.tbl,CERAtotalQUAN.black.adjusted.tbl),
  group_header = c("Dietary.LZ", "Dietary.LZ","Total L and Z","Total L and Z","Quartile of total L and Z","Quartile of total L and Z")) 

#CERAD: Trial 1 Score 
CERAtotal.White.adjusted.tbl <- svyglm(CERAD1 ~ Dietary.LZ+ Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment,design = subset(NHANES_design, Race == 'Non-Hispanic White')) %>%
  tbl_regression(exponentiate = TRUE, include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'White'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")

CERAtotal.black.adjusted.tbl <- svyglm(CERAD1 ~ Dietary.LZ+ Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment ,design = subset(NHANES_design, Race == 'Non-Hispanic Black')) %>%
  tbl_regression(exponentiate = TRUE, include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'Black'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
# Total L and Z
CERAtotalLZ.White.adjusted.tbl <- svyglm(CERAD1 ~ Total.LZ+ Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment,design = subset(NHANES_design, Race == 'Non-Hispanic White')) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'White'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
CERAtotalLZ.black.adjusted.tbl <- svyglm(CERAD1 ~ Total.LZ+ Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment ,design = subset(NHANES_design, Race == 'Non-Hispanic Black')) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'Black'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
#Quartile of total L and Z
CERAtotalQUAN.White.adjusted.tbl <- svyglm(CERAD1 ~ Total.LZ.quantile.var+ Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment,design = subset(NHANES_design, Race == 'Non-Hispanic White')) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'White'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
CERAtotalQUAN.black.adjusted.tbl <- svyglm(CERAD1 ~ Total.LZ.quantile.var+ Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment ,design = subset(NHANES_design, Race == 'Non-Hispanic Black')) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'Black'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
table5cerad1 <- tbl_stack(
  tbls = list(CERAtotal.White.adjusted.tbl, CERAtotal.black.adjusted.tbl,CERAtotalLZ.White.adjusted.tbl,CERAtotalLZ.black.adjusted.tbl,CERAtotalQUAN.White.adjusted.tbl,CERAtotalQUAN.black.adjusted.tbl),
  group_header = c("Dietary.LZ", "Dietary.LZ","Total L and Z","Total L and Z","Quartile of total L and Z","Quartile of total L and Z")) 


#CERAD: Trial 2 Score
CERAtotal.White.adjusted.tbl <- svyglm(CERAD2 ~ Dietary.LZ+ Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment,design = subset(NHANES_design, Race == 'Non-Hispanic White')) %>%
  tbl_regression(exponentiate = TRUE, include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'White'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")

CERAtotal.black.adjusted.tbl <- svyglm(CERAD2 ~ Dietary.LZ+ Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment ,design = subset(NHANES_design, Race == 'Non-Hispanic Black')) %>%
  tbl_regression(exponentiate = TRUE, include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'Black'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
# Total L and Z
CERAtotalLZ.White.adjusted.tbl <- svyglm(CERAD2 ~ Total.LZ+ Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment,design = subset(NHANES_design, Race == 'Non-Hispanic White')) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'White'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
CERAtotalLZ.black.adjusted.tbl <- svyglm(CERAD2 ~ Total.LZ+ Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment ,design = subset(NHANES_design, Race == 'Non-Hispanic Black')) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'Black'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
#Quartile of total L and Z
CERAtotalQUAN.White.adjusted.tbl <- svyglm(CERAD2 ~ Total.LZ.quantile.var+ Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment,design = subset(NHANES_design, Race == 'Non-Hispanic White')) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'White'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
CERAtotalQUAN.black.adjusted.tbl <- svyglm(CERAD2 ~ Total.LZ.quantile.var+ Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment ,design = subset(NHANES_design, Race == 'Non-Hispanic Black')) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'Black'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
table5cerad2 <- tbl_stack(
  tbls = list(CERAtotal.White.adjusted.tbl, CERAtotal.black.adjusted.tbl,CERAtotalLZ.White.adjusted.tbl,CERAtotalLZ.black.adjusted.tbl,CERAtotalQUAN.White.adjusted.tbl,CERAtotalQUAN.black.adjusted.tbl),
  group_header = c("Dietary.LZ", "Dietary.LZ","Total L and Z","Total L and Z","Quartile of total L and Z","Quartile of total L and Z")) 

#CERAD: Trial 3 Score
CERAtotal.White.adjusted.tbl <- svyglm(CERAD3 ~ Dietary.LZ+ Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment,design = subset(NHANES_design, Race == 'Non-Hispanic White')) %>%
  tbl_regression(exponentiate = TRUE, include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'White'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")

CERAtotal.black.adjusted.tbl <- svyglm(CERAD3 ~ Dietary.LZ+ Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment ,design = subset(NHANES_design, Race == 'Non-Hispanic Black')) %>%
  tbl_regression(exponentiate = TRUE, include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'Black'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
# Total L and Z
CERAtotalLZ.White.adjusted.tbl <- svyglm(CERAD3 ~ Total.LZ+ Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment,design = subset(NHANES_design, Race == 'Non-Hispanic White')) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'White'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
CERAtotalLZ.black.adjusted.tbl <- svyglm(CERAD3 ~ Total.LZ+ Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment ,design = subset(NHANES_design, Race == 'Non-Hispanic Black')) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'Black'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
#Quartile of total L and Z
CERAtotalQUAN.White.adjusted.tbl <- svyglm(CERAD3 ~ Total.LZ.quantile.var+ Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment,design = subset(NHANES_design, Race == 'Non-Hispanic White')) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'White'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
CERAtotalQUAN.black.adjusted.tbl <- svyglm(CERAD3 ~ Total.LZ.quantile.var+ Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment ,design = subset(NHANES_design, Race == 'Non-Hispanic Black')) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'Black'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
table5cerad3 <- tbl_stack(
  tbls = list(CERAtotal.White.adjusted.tbl, CERAtotal.black.adjusted.tbl,CERAtotalLZ.White.adjusted.tbl,CERAtotalLZ.black.adjusted.tbl,CERAtotalQUAN.White.adjusted.tbl,CERAtotalQUAN.black.adjusted.tbl),
  group_header = c("Dietary.LZ", "Dietary.LZ","Total L and Z","Total L and Z","Quartile of total L and Z","Quartile of total L and Z")) 

alllei5<-tbl_merge(tbls = list(table5sdr,table5cerad1,table5cerad2,table5cerad3),tab_spanner = c('CERAD: Score Delayed Recall','CERAD: Trial 1 Score','CERAD: Trial 2 Score','CERAD: Trial 3 Score'))
alllei5%>%as_flex_table()%>%flextable::save_as_html(path = "2.8table5.html")


#### tab5 结束####

#### tab6 开始####

# Dietary.LZ

tb5drwhite.age.adjusted.tbl <- svyglm(CERAD.total ~ Dietary.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, design = subset(NHANES_design, Race == 'Non-Hispanic White'))%>%
  tbl_regression(include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'White'),
                 pvalue_fun = function(x) style_pvalue(x, digits = 2)
  )%>%
  modify_header(
    std.error = "**SE**"
  )

tb5drwhite.age.adjusted.tbl

tb4dr.fully.adjusted.tbl <- svyglm(CERAD.total ~ Dietary.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                   design = NHANES_design) %>%
  tbl_regression( include = c(Dietary.LZ),
                  label = list(Dietary.LZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
# Supplement.LZ
tb4supp.age.adjusted.tbl <- svyglm(CERAD.total ~ Supplement.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Supplement.LZ),
                 label = list(Supplement.LZ ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
tb4supp.fully.adjusted.tbl <- svyglm(CERAD.total ~ Supplement.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                     design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Supplement.LZ),
                 label = list(Supplement.LZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
# Total L.Z gcr
tb4Total.age.adjusted.tbl <- svyglm(CERAD.total ~ Total.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
tb4Total.fully.adjusted.tbl <- svyglm(CERAD.total ~ Total.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                      design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
# Total L.Z Quartile gcr
tb4Totallz.age.Quartile.adjusted.tbl <- svyglm(CERAD.total ~ Total.LZ.quantile.var + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
tb4Totallz.fully.adjusted.tbl <- svyglm(CERAD.total ~ Total.LZ.quantile.var + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                        design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'Fully-adjusted'))
# regression 的结果一般采用纵向拼接（stack）而不是横向拼接（merger）
tb4CERAD.total<-tbl_stack(
  tbls = list(tb4dr.age.adjusted.tbl, tb4dr.fully.adjusted.tbl,
              tb4supp.age.adjusted.tbl, tb4supp.fully.adjusted.tbl,tb4Total.age.adjusted.tbl,tb4Total.fully.adjusted.tbl,tb4Totallz.age.Quartile.adjusted.tbl,tb4Totallz.fully.adjusted.tbl),
  group_header = c("Dietary.LZ", "Dietary.LZ", "Supplement.LZ", "Supplement.LZ","Total.LZ", "Total.LZ","Quartile of total L and Z","Quartile of total L and Z")) 

# Animal Fluency score
Animaldr.age.adjusted.tbl <- svyglm(Animal.Fluency ~ Dietary.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
Animaldr.fully.adjusted.tbl <- svyglm(Animal.Fluency ~ Dietary.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                      design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
# Supplement.LZ
Animalsupp.age.adjusted.tbl <- svyglm(Animal.Fluency ~ Supplement.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Supplement.LZ),
                 label = list(Supplement.LZ ~ 'Age-adjusted'))
Animalsupp.fully.adjusted.tbl <- svyglm(Animal.Fluency ~ Supplement.LZ + Age + Sex + BMI + Alq.group + 
                                          Smoke.group + PIR + Education.attainment, 
                                        design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Supplement.LZ),
                 label = list(Supplement.LZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )

# Total L.Z gcr
AnimalTotal.age.adjusted.tbl <- svyglm(Animal.Fluency ~ Total.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
AnimalTotal.fully.adjusted.tbl <- svyglm(Animal.Fluency ~ Total.LZ + Age + Sex + BMI + Alq.group + 
                                           Smoke.group + PIR + Education.attainment, 
                                         design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
# Total L.Z Quartile gcr
AnimalTotallz1.age.Quartile.adjusted.tbl <- svyglm(Animal.Fluency ~ Total.LZ.quantile.var + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
AnimalTotallz1.fully.adjusted.tbl <- svyglm(Animal.Fluency ~ Total.LZ.quantile.var + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                            design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
# regression 的结果一般采用纵向拼接（stack）而不是横向拼接（merger）
Animal.Fluencytab<-tbl_stack(
  tbls = list(Animaldr.age.adjusted.tbl, Animaldr.fully.adjusted.tbl,
              Animalsupp.age.adjusted.tbl, Animalsupp.fully.adjusted.tbl,AnimalTotal.age.adjusted.tbl,AnimalTotal.fully.adjusted.tbl,AnimalTotallz1.age.Quartile.adjusted.tbl,AnimalTotallz1.fully.adjusted.tbl),
  group_header = c("Dietary.LZ", "Dietary.LZ", "Supplement.LZ", "Supplement.LZ","Total.LZ", "Total.LZ","Quartile of total L and Z","Quartile of total L and Z")) 


# Digit Symbol Score
#Dietary L and Z
dss.age.adjusted.tbl <- svyglm(DSST ~ Dietary.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )

dss.fully.adjusted.tbl <- svyglm(DSST ~ Dietary.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                 design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )

# Supplement.LZ
dsssupp2.age.adjusted.tbl <- svyglm(DSST ~ Supplement.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Supplement.LZ),
                 label = list(Supplement.LZ ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
dsssupp2.fully.adjusted.tbl <- svyglm(DSST ~ Supplement.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                      design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Supplement.LZ),
                 label = list(Supplement.LZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
# Total L.Z gcr
dssTotal2.age.adjusted.tbl <- svyglm(DSST ~ Total.LZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
dssTotal2.fully.adjusted.tbl <- svyglm(DSST ~ Total.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                       design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ),
                 label = list(Total.LZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )

# Total L.Z Quartile gcr
dssTotallz2.age.Quartile.adjusted.tbl <- svyglm(DSST ~ Total.LZ.quantile.var + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )

dssTotallz2.fully.adjusted.tbl <- svyglm(DSST ~ Total.LZ.quantile.var + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
                                         design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(Total.LZ.quantile.var),
                 label = list(Total.LZ.quantile.var ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
# regression 的结果一般采用纵向拼接（stack）而不是横向拼接（merger）
DSSTcol<-tbl_stack(
  tbls = list(dss.age.adjusted.tbl, dss.fully.adjusted.tbl,
              dsssupp2.age.adjusted.tbl, dsssupp2.fully.adjusted.tbl,dssTotal2.age.adjusted.tbl,dssTotal2.fully.adjusted.tbl,dssTotallz2.age.Quartile.adjusted.tbl,dssTotallz2.fully.adjusted.tbl),
  group_header = c("Dietary.LZ", "Dietary.LZ", "Supplement.LZ", "Supplement.LZ","Total.LZ", "Total.LZ","Quartile of total L and Z","Quartile of total L and Z")) 
# regression 的结果一般采用纵向拼接（stack）而不是横向拼接（merger）

alllei4<-tbl_merge(tbls = list(tb4CERAD.total,Animal.Fluencytab,DSSTcol),tab_spanner = c('CERAD: Total score','Animal Fluency score','Digit Symbol Score'))
alllei4%>%as_flex_table()%>%flextable::save_as_html(path = "2.8table4.html")



#### tab6 结束####


show_header_names(alllei)
alllei%>%  modify_column_hide(p.value)%>% # 隐藏p值
  modify_header(estimate ~ 'exp(Beta) (95% CI)')%>% # 编辑表头
  modify_column_merge(pattern = "{estimate} ({ci})")  # 合并列值

as_flex_table(alllei) %>%
  flextable::save_as_html(path = 'aaatbl_stack_result.html')


#### 4. Subgroup ####
# # 以性别为例，区分男性以及女性
# # 男性
male.supp.age.adjusted.tbl <- svyglm(CERAD.delay.recall ~ Supplement.LZ + Age,
                                design = subset(NHANES_design, Sex == 'male')) %>% # 注意，在这里用 subset 进行选取
  tbl_regression(exponentiate = TRUE, include = c(Supplement.LZ),
                 label = list(Supplement.LZ ~ 'Age-adjusted'))
# 
# # 女性
# female.supp.age.adjusted.tbl <- svyglm(CERAD.delay.recall ~ Supplement.LZ + Age, 
#                                      design = subset(NHANES_design, Sex == 'female')) %>% # 注意，在这里用 subset 进行选取
#   tbl_regression(exponentiate = TRUE, include = c(Supplement.LZ),
#                  label = list(Supplement.LZ ~ 'Age-adjusted'))
# 
# # 拼接
# tbl_stack(
#   tbls = list(male.supp.age.adjusted.tbl, female.supp.age.adjusted.tbl),
#   group_header = c("Male", "Female")) 
# 
# 
# #### 5. 结果输出为 word ####
# tbl_stack(
#   tbls = list(dr.age.adjusted.tbl, dr.fully.adjusted.tbl,
#               supp.age.adjusted.tbl, supp.fully.adjusted.tbl),
#   group_header = c("Dietary.LZ", "Dietary.LZ", "Supplement.LZ", "Supplement.LZ")) %>%
#   as_flex_table() %>%
#   flextable::save_as_docx(path = 'tbl_stack_result.docx')


  